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ABSTRACT 

We present the Multifrequency Snapshot Sky Survey (MSSS), the first northern-sky LOFAR imaging survey. In this introductory paper, we first 
describe in detail the motivation and design of the survey. Compared to previous radio surveys, MSSS is exceptional due to its intrinsic multi¬ 
frequency nature providing information about the spectral properties of the detected sources over more than two octaves (from 30 to 160 MHz). 
The broadband frequency coverage, together with the fast survey speed generated by LOFAR’s multibeaming capabilities, make MSSS the first 
survey of the sort anticipated to be carried out with the forthcoming Square Kilometre Array (SKA). Two of the sixteen frequency bands included 
in the survey were chosen to exactly overlap the frequency coverage of large-area Very Large Array (VLA) and Giant Metrewave Radio Telescope 
(GMRT) surveys at 74 MHz and 151 MHz respectively. The survey performance is illustrated within the “MSSS Verification Field” (MVF), a re¬ 
gion of 100 square degrees centered at (a, S)j 2 ooo = (15^, 69°). The MSSS results from the MVF are compared with previous radio survey catalogs. 
We assess the flux and astrometric uncertainties in the catalog, as well as the completeness and reliability considering our source finding strategy. 
We determine the 90% completeness levels within the MVF to be 100 mJy at 135 MHz with 108" resolution, and 550 mJy at 50 MHz with 166" 
resolution. Images and catalogs for the full survey, expected to contain 150000-200000 sources, will be released to a public web server. We 
outline the plans for the ongoing production of the final survey products, and the ultimate public release of images and source catalogs. 

Key words. Surveys — Radio continuum: general 


1. Background 

All-sky continuum surveys are a key application of radio tele¬ 
scopes. They provide a view of galaxies across the Universe free 
from the biasing effect of extinction, enable searches for rare 
sources, and provide a pathway for the discovery of new phe¬ 
nomena. Several large-area surveys have been performed with 
existing radio telescopes at a number of frequencies. Among 
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these, many of the earliest were performed at low frequen- 


cies (y < 350 MHz) (see e.g. [Jauncey[[1975 

). The new Low 

Frequency Array (LOFAR; [van Haarlem et a 

[.[[2013J operates 


at frequencies between 10 and 240 MHz, and was constructed 
with one primary aim being to perform groundbreaking imaging 
surveys of the northern sky ( [Rottgering et al.||2011| ). Currently, 
the most extensive low frequency catalogs at or near LOFAR 
frequencies are the Eighth Cambridge Survey of Radio Sources 
(8C; |Rees||199^ catalog revised by [Hales et ~ar]|1995|), the 
VLA Low-Frequency Sky Survey (VLSS; [Cohen et al.||2007 ' 
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and the revised catalog VLSSr; Lane et al.||20] 

[2|, the Seventh 

Cambridge Survey of Radio Sources (7C; [Ha 

es et al.||2007). 

the Westerbork Northern Sky Survey (WENSS; 

Rengelink et a . 

1997]), and most recently the Murchison Widefie 

id Array (MWA; 


Tingay et al 


Walker et al. 


|2Q13| ) Commissioning Survey (MWACS; |Hurley-| 


2015| ). Other ongoing low-frequency radio surveys 


include the TIFR GMRT Sky Survey (TGSSQ) that is using the 
Giant Metrewave Radio Telescope (GMRT) and has already re¬ 
leased survey data products to the commu nity, and the Galac tic 
and Extragalactic MWA Survey (GLEAM; |Wayth et al.|2015| ). 

A key application of all-sky radio surveys is the comparison 
of source properties at the wide range of frequencies at which 
they are detected. This provides crucial information from which 
the physical properties of these sources can be derived. 

To date, no wide-area radio surveys have been performed 
with large fractional bandwidth (i.e., 2:1 or more). This situa¬ 
tion is bound to change in the era of the Square Kilometre Array 
(SKA; |Carilli & Rawlings]|2004| ). The SKA and its pathfinder 
and precursor projects plan wide-area surveys of radio con¬ 
tinuum sources, with large fractional bandwidth. This opens 
the door for the spectral study of sources detected within the 
survey, using only the survey data itself. LOEAR is a key 
SKA pathfinder telescope ( |van Haarlem et al!]|2013| ) in the 10 
and 240 MHz frequency range. The array is centered in the 
Netherlands with current outlying stations in Germany, Erance, 
the United Kingdom, and Sweden. LOEAR is built up from thou¬ 
sands of dipoles clustered in groups called stations. The signals 
from the dipoles making up each station are digitally combined 
to steer the beam in one or more directions of interest. Stations 
are combined in a software correlator located in Groningen, a 
city in the north of the Netherlands. 


One of the key applications of LOEAR is wide-field imag¬ 
ing. In this paper we introduce a new radio survey performed 
with LOEAR, the Multifrequency Snapshot Sky Survey (MSSS), 
that has been driven forward as a commissioning project for the 
telescope. The MSSS project serves as a testbed for operations - 
particularly large-scale imaging projects - and enables straight¬ 
forward processing of later observations. 


The motivations for performing MSSS and the design of the 
survey are described in detail in §|^ The calibration and imag¬ 
ing strategies are presented §[3) and the resulting standard data 
products are described in §|4pWe illustrate the performance of 
the survey through a detailed analysis of the “MSSS Verification 
Eield” (MVE) in §[^ Several avenues for scientific exploitation 
of MSSS are outlined in §[^ and we conclude the paper in §[ 7 ] 


2. Context and survey design 

Imaging applications with the LOEAR telescope will require au¬ 
tomated processing. The calibration step in particular needs to 
be largely unattended, with a major implication that a priori sky 
models are required at arbitrary locations on the sky. A number 
of northern sky radio surveys are available in the literature, but 
do not cover the proper frequency range at the resolution needed 
to reliably initiate the calibration routines. Moreover, a coher¬ 
ent commissioning project was required to focus development 
activities and produce a generic automated processing pipeline, 
while simultaneously exercising the end-to-end telescope oper- 


^ http://tgss.ncra.tifr.res.in/150MHz/tgss.html 


ations. These goals led to the initiation of the “Multifrequency 
Snapshot Sky Survey” (MSSS|^. 

In the original MSSS plans (early 2008), it was anticipated 
that the survey would be performed using 13 core stations 
(CS), 7 remote stations (RS), and 3 international stations (IS). 
Ultimately, the array construction proceeded rapidly, and MSSS 
has been performe d wit h the full complement of stations (ex¬ 
cept in HBA; see §^). This leads to a significantly different 
array than originally envisioned, both in terms of sensitivity and 
uv coverage. A complete overview of the LOEAR system is pro¬ 
vided by van Haarlem et al. ( 2013 ). The telescope layout, as well 
as processing software limitations (now substantially reduced), 
led to plans for a low-resolution survey initially, with aspirations 
for a higher resolution survey in future. This paper represents 
the initial low-resolution effort, and we present prospects for our 
plans for higher resolution data products in § |6.1[ 

The MSSS survey effort needs to provide images and cata¬ 
logs with sufficient angular resolution to reliably initiate the self¬ 
calibration cycle in the imaging pipeline. At the same time, the 
frequency coverage needs to be sufficient to ensure that spectral 
variations within the model are accounted for. These require¬ 
ments were balanced with the need for a relatively rapid survey, 
taking on the order of weeks of telescope time to perform. 

With all of these considerations in mind, we designed a 
two-component observational strategy. The Low Band Antenna 
(LBA) component covers the 30-75 MHz range, and the High 
Band Antenna (HBA) component covers the 119-158 MHz 
range. The exact frequencies were chosen to evenly cover the 
LB A range, and to avoid major radio frequency interference 
(REI) in the HBA range (see Table [^. The number of frequency 
bands (eight 2 MHz bands in each of the LBA and HBA com¬ 
ponents) were chosen to allow multiplexing the sky coverage. 
In early survey test observations, the “16-bit” correlator mode 
allowed three fields to be observed simultaneously, each with 
16 MHz bandwidth. Near the end of 2012 (when most of the 
LBA test observations were complete, and test HBA observa¬ 
tions were beginning), the “8-bit” correlator mode became avail¬ 
able, doubling the number of fields that can be simultaneously 
observed (each with 16 MHz bandwidth) to six. All observations 
provide data in all four Stokes parameters. The key parameters 
for the two frequency components of the survey are summarized 
in Table and the setup of these are described in turn below. 


2.1. Setup of MSSS-LBA 

The LBA component of MSSS is observed using the 
LBAJNNER configuration. In this mode, the digital station 
beams are formed using signals from the inner 48 dipoles of 
each 96-dipole station. The resulting compact station, with di¬ 
ameter T) = 32.25 m, provides a large field of view. International 
stations are included in all MSSS-LBA observations in addi¬ 
tion to the Dutch stations, and these come with the full com¬ 
plement of 96 dipoles. The LBA survey pointings were designed 
using a nominal primary beam half-power beam width (HPBW) 
at 60 MHz of 1U55 (from HPBW = using = 1.3 

and A = 5 m). We now know that the appro priate value of 
at = 1.10 + 0.02 ( [van Haarlem et al^|2013| ), so the station 
beams are 15% smaller than initially anticipated, and a some¬ 
what larger variation in image noise across the survey can there¬ 
fore be expected. This will be far less evident at MSSS-LBA 


^ The original name of the survey was the “Million Source Shallow 
Survey”. Under current projections (see §[^, we expect to catalog well 
over 10^ sources, but probably not 10^. 
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frequencies lower than 60 MHz, where the beam sizes are much 
larger. Survey fields are laid out on strips of constant declination, 
and spaced equally on each strip by Aa < HPBW/2. The decli¬ 
nation strips are themselves separated by AS = HPBW/2. This 
results in a total of 660 MSSS-LBA fields. 

Initial tests, performed before the station digital beam form¬ 
ing was properly calibrated, indicated that interleaving calibrator 
observations between target observations would not lead to a suf¬ 
ficiently stable amplitude scale (see §[^for details of the calibra¬ 
tion strategy). Thus, we adopted an observing mode wherein one 
of the primary calibrators is always observed, using the identi¬ 
cal 16 MHz bandwidth, in parallel with two simultaneous target 
field observations (or five target fields, when the 8-bit mode is 
used). 

The amount of observing time used per field is based on the 
goal of obtaining image noise at the 10 mJy beam“^ level, in each 
of the eight 2 MHz bands. Using early projections of the image 
sensitivity expected from LOFAR, an estimate of 90 minutes per 
field was obtained. To improve uv coverage (see below) this was 
split up into 9 snapshots. The start times of the individual snap¬ 
shots are spaced by 1 h for northern target fields (S > 30°), 
and by 45 min for target fields closer to the celestial equator 
(S < 30°). The central (fifth) snapshot is observed near transit 
(HA 0 h). For ease of scheduling, the final survey observing 
time per field is 9 x 11 =99 min, yielding initially estimated the¬ 
oretical thermal noise values between about 6-20 mJybeam“^ 
depending on band. 

From observations of calibrator sources, we now have em¬ 
pirical estimates of the LOFAR station Sys tem Equivalent Flux 
Densit ies (SEFDs); these are provided by |van Haarlem et al^ 
( |2013| ). Based on those numbers we can see that an observing 
time of 99 min yields an expected thermal noise between about 
5-10 mJy beam“^. It should be noted however that the rudimen¬ 
tary calibration procedures that are implemented in the MSSS 
pipeline limit our actual sensitivity to an image noise that is typ¬ 
ically a factor of a few higher than the thermal noise estimate. 
Moreover, classical confusion noise would likely limit images 
produced at the limited angular resolution targeted for default 
MSSS imaging (out to uv distance of 2-3 kT), see Fig. We 
calculate the expected confusion noise in two ways. The first is 
based on extrapolation f rom VLSS B- configuration estimates of 
the confusion noise (see |Cohen|2004| ) : 

/ 0 / y 

crconf,VLSs=29(-J //Jybeam-, (1) 

where 0 is the synthesized beamsize, v is the observing fre¬ 
quency, and we have extrapolated the VLSS estimate usi ng a typ- 
ical spectral index of -0.7. The second estimate is from [Condon 
et al.| ( |2Qf^ , who used deep VLA C-configuration observations 
at S-band (2-4 GHz) to derive 

/ 0 / y \-0-7 

O-conf,Condon = L2 — j ( 3 q 2 QHz / 

We consider the numbers predicted for MSSS fairly reliable 
since the two estimates provide very similar values (on average 
only 6.4% different across the full MSSS frequency range) de¬ 
spite being based on data from very different observing frequen¬ 
cies and angular resolutions. Still, the sensitivity obtained in the 
MSSS data presented in §[^ suggests that the confusion limit is 
somewhat lower than predicted, at least in the HBA. 

The uv coverage obtained by splitting the observations of 
LB A fields into 9 snapshots is illustrated in Fig.The figure 
was created using the observations described in §[^ 
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Fig. 1. Theoretical thermal and confusion noise per band for 
MSSS observations. The valu es take into account empir ically- 
determined SEFD values from [van Haarlem et al.| ( |20l'3] ), and a 
uv distance cutoff of 3 kT. 


Including overheads, the total amount of observing time re¬ 
quired to complete MSSS-LBA (assuming 8-bit mode) is 297 h. 
Approximately 130 TB of raw visibility data is collected by ob¬ 
serving all 660 pointings and associated calibrator scans. 

During the LBA portion of the survey, a single subband 
(bandwidth 195 kHz) centered at 60 MHz is always placed on 
the north celestal pole (NCP) as part of a transient and variabil¬ 
ity monitoring campaign coordinated by the LOFAR Transients 
Key Project. We return to this feature of the survey setup in §[^ 


2.2. Setup of MSSS-HBA 


The HBA component of MSSS is observed using the 
HBA_DUAL_INNER station configuration. In this mode, both 
24-tile substations are utilized separately for each of the core 
stations. The outer 24 tiles of the 48-tile remote stations are dis¬ 
abled so that the field of view of the stations are all identical (at a 


small cost i n sensitivity on the longer baselines; see van Haarlem 
et al.||2013 ). The survey pointings were designed using a nom¬ 
inal beam size at 150MHz of HPBW = 4?84 (using = 1.3, 
T = 2 m and D = 30.75 m). As in the LBA portion of the sur¬ 
vey, our understanding of the station beam sizes has now been 
empirically determined to be smaller (a^ = 1.02 + 0.01 for HBA 
core stations; [van Haarlem et al.|2013| ). International stations are 
not included in the HBA portion of the survey, because at the 
time of observations the system had restrictions on the number 
of correlator inputs, requiring the loss of core stations when in¬ 
ternational stations were included. Since MSSS is primarily a 
low angular resolution survey, we kept the full complement of 
core stations in the survey observations. A separate and com¬ 
plementary survey including the international stations has been 
conducted to search a representativ e portion of the LOFA R HBA 
sky for compact (< 1") calibrators ( Moldon et al.|2Q15 ). 

With the HBA observations, a simultaneous calibrator strat¬ 
egy as implemented in MSSS-LBA is not required, since the sta¬ 
bility is much better. Moreover, such a strategy would not be pos¬ 
sible, because the analog beamformer in each of the 16-dipole 
tiles reduces the field of view to approximately HPBWtiie = 
A/Dhiq = 22.9 deg (where Ane = 5 m is the tile size) at 150 MHz. 
Thus, most field observations would not be able to include a par¬ 
allel primary calibrator observation with sufficient sensitivity. 
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Table 1. LOFAR MSSS configuration 


Parameter 

MSSS-LBA 

MSSS-HBA 

Station configuration 

Band central frequencies (MHz) 

LBAJNNER 

HBA_DUAL_INNER 

Band 0 

31 

120 

Band 1 

37 

125 

Band 2 

43 

129 

Band 3 

49 

135 

Band 4 

54 

143 

Band 5 

60 

147 

Band 6 

66 

151 

Band 7 

74 

157 

Calibrator observation type 

Simultaneous 

Serial 

Calibrator observation time (min) 

11 

1 

Observing time per snapshot (min) 

11 

7 

Number of snapshots 

9 

2 

Snapshot gap (h) 

1 (S > 30°) 

0.75 (S < 30°) 

4 

Number of fields 

660 

3616 

Total survey time with overheads (h) 

297 

201 
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Fig. 2. The uv coverage for the MSSS-LBA field L227+69, core only (left) and all Dutch stations (right). 


Calibrator observations are therefore performed between field 
observations, using the same bandwidth as the field observa¬ 
tions. 


Similar sensitivity considerations as for the LBA component 
of MSSS led to a required observing time per HBA field of ap¬ 
proximately 15 min. Snapshots are also used for HBA obser¬ 
vations to somewhat improve uv coverage. Two snapshots are 
used, with the start times of the snapshots separated by 4 h 
(bracketing transit, such that the hour angles of the snapshots 
are HA +2 h). For ease of scheduling we adopted 2x7 min 
integrations per field. This leads to an estimated thermal noise of 
about 1 mJy beam“^ per band, using the empirical SEFD values 
given by [van Haarlem et al. ( [2013| ). As with the LBA sensitiv¬ 
ity, the actual noise level in HBA images is a factor of a few 
higher than the thermal noise estimate, and confusion is likely 
the true limiting factor in images with limited angular resolu¬ 
tion, see Fig.[2 


The uv coverage of the two-snapshot HBA observing strat¬ 
egy is shown in Fig.[^ The survey grid was designed in the same 
way as the LBA grid, using HPBW=4?84. The HBA component 
of MSSS is made up of 3616 fields. Including overheads, the to¬ 
tal amount of observing time required to complete MSSS-HBA 
(assuming 8-bit mode) was 201 h. Approximately 470 TB of 
raw visibility data is collected by observing all 3616 pointings 
and associated calibrator scans. 

2.3. Survey fields 

The layout of the MSSS survey fields was determined in such 
a way as to provide nearly uniform coverage at the optimized 
frequencies 60 MHz (MSSS-LBA) and 150 MHz (MSSS-HBA), 
as described in § |2.1| and [2!^ The coordinates of the center of the 
survey fields are shown on an Aitoff projection in Fig.|^ 

The MSSS fields planned and observed so far have a lower 
declination limit of ^ = 0°. At a later date, the lower declination 
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Fig. 3. The uv coverage for the MSSS-HBA field H229+70, core only (left) and all Dutch stations (right). 



75 ° 




Fig. 4. MSSS fields for the LB A (top) and HBA (bottom) por¬ 
tions of the survey, presented in equatorial coordinates on an 
Aitoff projection. The LB A survey consists of 660 fields, while 
the HBA survey consists of 3616 fields. 

limit may be extended farther to the south, extending the MSSS 
sky coverage. 

2.4. Survey mosaics 

The final presentation of the survey images (and derivation of 
the resulting catalog) will be based on 10° x 10° mosaics gener¬ 
ated from the individual HBA and LB A fields. The mosaic grid 


is common for both bands in order to facilitate multi-frequency 
fiux comparison of sources detected in the survey. The survey 
contains a total of 214 mosaic fields, including one NCP mosaic 
centered at ^ = 90°. 


3. Calibration and imaging strategy 

The MSSS commissioning project has driven the development 
of the first version of the standard imaging pipeline (SIP; |HeaTd 
et al.|2010| ), which can be scheduled by the control system to run 
automatically on the central processing cluster upon completion 
of the individual observations. The SIP embodies the calibration 
strategy described in this section, and is now being upgraded to 
allow improved image quality at higher angular resolution. 

We illustrate the processing chain followed by a single ob¬ 
servation in Fig. Fig. shows a schematic view of the over¬ 
all processing pipeline, which combines multiple snapshots and 
calibration observations to create final combined images and 
feeds source catalogs into the LOFAR Global Sky Model (GSM) 
database. In these diagrams, cylinders indicate data products, 
and blocks indicate programs and pipeline segments. The full 
pipeline is broken up into three main components. The first is 
the Calibrator Pipeline (CP), which applies some pre-processing 
steps that are described below to the calibrator scans, and per¬ 
forms primary calibration using a known and well-understood 
calibrator model. The Target Pre-processing Pipeline (TPP) per¬ 
forms the same pre-processing steps as in CP, and subsequently 
uses the station gains derived in CP to apply the primary cali¬ 
bration to the individual field snapshots. The resulting calibrated 
snapshots are stored until all snapshots of a particular field are 
completed, after which the Target Imaging Pipeline (TIP) be¬ 
gins. This final stage is the heart of the pipeline, and consists of 
imaging the field and optionally running a self-calibration cycle 
which is still undergoing further development. We now proceed 
to detail the various pieces of the pipeline. 
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Fig. 5. Sketch of the MSSS calibration and imaging pipeline. See the text for a description. 



Fig. 6. The scheme for calibrating, processing, and combining 
individual snapshots of an MSSS field. Each field is observed 
in N snapshots, and each snapshot observation simultaneously 
observes M beam directions. In the case of MSSS-LBA, N = 9, 
and M is either 2 or 5 depending on whether the observations 
are done in 16-bit or 8-bit mode. For MSSS-HBA, N = 2 and 
M = 3 or 6. CP stands for Calibrator Pipeline; TPP for Target 
Pre-processing Pipeline; TIP for Target Imaging Pipeline; and 
GSM for Global Sky Model. See Fig. for an overview of how 
these pieces fit together in more detail, and for the steps that 
make up each segment of the full pipeline. 


3.1. Pre-processing steps 


Following the fiagging step, the demixing technique ( |van| 
|der Tol et al.||2007] ) is applied in order to remove far off-axis, 
bright sources (primarily the so-called “A-team”: Cygnus A, 
Cassiopeia A, Virgo A, and Taurus A) from the visibility data. 
The automatic pipeline calculates the distance to the A-team 
sources from the target (and calibrator) fields. Sources within 
distance ranges determined empirically from commissioning ex¬ 
perience are selected for demixing. 

Compression of the data to a manageable volume is done au¬ 
tomatically following the demixing step. Typically, in both LB A 
and HBA the data volume is reduced by a factor of 10 after av¬ 
eraging and calibration have been performed. The compression 
factors were selected to minimize bandwidth and time smearing 
effects, as well as retain sufficient time resolution to allow cap¬ 
turing time variable ionospheric effects. In this way we are able 
to restrict our estimates of position-dependent survey sensitiv¬ 
ity to only consider the combination of station beam pattern and 
survey pointing grid. 


3.1.1. Bandwidth smearing 


The effect of finite bandwidth is to partially decorrelate the sig¬ 
nal and leads to a radial smearing of sources far from the phase 
tracking center. Assuming a square bandpass and a synthesized 
beam with a Gaussian profile, the magnitude of the reduction 


in peak flux de nsity can be approximated as given by Bridle & 
|Schwab| ( |l999| ): 


Flagging for RFI is performed in a standard and automatic fash¬ 
ion using the AOFlagger ( [Qffringa et al.|201()| ). This program has 
been shown to provide excellent RFI excision with minimal false 
positives. Details of the typical performance of the implemented 
algorithm on LOFA R data, along with repr esentative RFI statis¬ 
tics, are presented by Qffringa et al.| ( |2013] ). 


I 

To 


SVc f 

2Vta2>-Av" 


OVc 


( 3 ) 


where 6 is the synthesized beam size (FWHM), is the central 
frequency of the observation, r is the angular distance from the 
phase center, and Av is the bandwidth. 
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For LOFAR, the individual subband width (using the stan¬ 
dard 200 MHz clock) is 195.3125 kHz and in MSSS observa¬ 
tions is divided into 64 channels. For a characteristic angular 
resolution of 2', we have calculated smearing factors for each 
of the MSSS band frequencies, and a field r adiu s co rresp ond- 
ing to half of the station beam HPBW (see § |2.1| and [Z^ . The 
resulting bandwidth smearing curves have been calculated for 
different frequency averaging parameters and are displayed in 
Fig-IZl 

In the case of the LBA survey, at least 8 channels per sub¬ 
band must be retained in order to keep the effect of smearing to 
<1-2%. For the HBA portion, the frequency averaging does not 
play such a crucial role. We thus retain 8 channels per subband 
in LBA for a maximum allowable bandwidth smearing factor at 
the lowest frequencies. In HBA, we retain 4 channels per sub¬ 
band to allow reprocessing at higher angular resolution without 
unacceptable smearing losses. Note that in the LBA, reprocess¬ 
ing the data in order to image at higher resolution will require 
either redoing the pre-processing steps with a lower frequency 
averaging factor, or acceptance of a substantial smearing factor 
at the lowest frequencies. 


3.1.2. Temporal resolution 

Another effect that needs to be accounted for is tim e-average 
smearing. Agai n referring to the expressions given by Bridle & 
Schwab] (1999| ), we assess the impact on MSSS data using: 


- = 1 - 1.22 X 10"^ 
fo 


(^r 


(4) 


where Ta is the averaging time. The estimated time smearing fac¬ 
tors are shown in Figure |7] 

The effect of time smearing is considerably less than the im¬ 
pact of bandwidth smearing. For small time averaging intervals, 
the effect is negligible. The smearing is also less important than 
the need to retain high time resolution for recovery of iono¬ 
spheric phase disturbances during the calibration process. We 
average the time steps to 2 seconds in order to recover high- 
quality station gain phases. This is illustrated in §[^ 


3.2. Primary and secondary calibration 

The primary calibrators for MSSS are listed in Table and are 
based on the calibration model presented by [Scaife & Heald 
( |2012| ). The tabulated spectral coefficients correspond to that for¬ 
mulation, namely the (n > 1) factors in 


log5'y(y) = logAo + Al logy -i- A2 log^ y -1-.. 


(5) 


These sources are mostly chosen to be compact on Dutch base¬ 
lines (approximately 80 km or less), bright enough to give suf¬ 
ficient signal-to-noise per visibility, and have well-characterized 
radio spectra. The exception is Cygnus A, which is used as a 
primary calibrator in the LBA portion of MSSS, despite being a 
very complicated extended source. It does however have a very 
well determined source model, based on extensive commission¬ 
ing work (summarized by McKean et al.|20lT| McKean|2015 ). 

Observational verification of these primary calibrator 
sources started early in the MSSS test program, and revealed 
that while the brightest sources in the low band (3C196, 3C295, 
and CygA) were suitable primary calibrators, the others were too 
weak and/or too close to A-team sources to provide stable gain 
amplitude solutions. These fainter sources are still useful for the 
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Fig. 7. Bandwidth (top) and time (bottom) smearing factors rel¬ 
evant to LOFAR MSSS observations, expressed as percentages 
and where unity is equivalent to no smearing. These values are 
calculated for 2' angular resolution and for sources at the half¬ 
power point of the station beam. In the top panel, points are plot¬ 
ted for averaging single subbands to 4 (stars), 8 (diamonds), and 
16 (squares) channels. In the bottom panel, the smearing is cal¬ 
culated for visibility time averaging intervals of 1 (squares), 2 
(diamonds), and 4 (stars) seconds. 


HBA portion of the survey but are not utilized in the LBA por¬ 
tion. 

Primary (flux) calibration is handled with two different 
strategies in the LBA and HBA parts of the survey. In the LBA 
part, we take advantage of the fact that the individual dipoles 
are sensitive to emission from the entire visible sky, and there is 
no analog beamformer limiting the field of view. This allows us 
to observe with a simultaneous calibrator beam. The calibrator 
at the highest elevation angle is used, regardless of the distance 
between calibrator and target fields. The calibrator beam uses 
the exact same frequency coverage as the target fields and runs 
for the full length (11 minutes) of the target snapshots. This en¬ 
sures sufficient time samples to obtain a stable gain amplitude 
solution. The CP obtains, and then exports, the median gain am¬ 
plitude solution per station, per snapshot (along with the corre¬ 
sponding gain phase that is not used) from the calibrator beam, 
for application to the flagged and demixed target data in the TPP. 
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Table 2. Primary MSSS calibrator sources 


Source ID 

RA (J2000.0) 

Dec (J2000.0) 

5 150 MHz (Jy) 

Spectral coefficients 

Morphology 

LBA calibrator 

3C48 

01'’37"’4E3 

+33°09'35" 

64.768 

(-0.387,-0.420,0.181) 

point 

No 

3C147 

05M2"’3ff.l 

+49°51'07" 

66.738 

(-0.022,-1.012,0.549) 

point 

No 

3C196 

08'’13"’36!0 

+48n3'03" 

83.084 

(-0.699,-0.110) 

double 

Yes 

3C286 

13'’31"’08!3 

+30°30'33" 

27.477 

(-0.158,0.032,-0.180) 

point 

No 

3C295 

14hlim2o?5 

+52n2'10" 

97.763 

(-0.582,-0.298,0.583,-0.363) 

double 

Yes 

3C380 

18'’29”’3E8 

+48M4'46" 

77.352 

(-0.767) 

point+diffuse 

No 

CygA 

19'’59”’28!3 

+40°44'02" 

10690.0 

(-0.670,-0.240,0.021) 

FRII 

Yes 


On the other hand, a different strategy is used for the HBA ob¬ 
servations. The analog HBA tile beam limits the field of view 
to typically 23 degrees HPBW at 150 MHz. This means that a 
bright, compact calibrator is not always available within the field 
of view near the targets. Therefore, the calibrator is observed 
alone (before the target snapshot) for 1 minute. Because the in¬ 
stantaneous sensitivity of the HBA system is much higher than 
that of the LBA system, stable gain amplitudes are obtained with 
a much shorter observing interval. As with LBA, the HBA cal¬ 
ibrator beam covers the same frequencies as the target beams, 
and the CP exports the median gain amplitude per station. These 
are subsequently applied to the target snapshot data in the TPP. 

Following the application of the primary flux calibration, the 
station phases remain uncalibrated (in the direction of the target 
field). The phase calibration takes place in the Target Imaging 


Pipeline, described in § 3.4 


3.3. Automatic data quality filtering 

Much of the MSSS-LBA data were obtained early in LOFAR’s 
lifetime, and this meant that much of the data were taken in un¬ 
stable conditions. It turned out that typical observations included 
one or more bad stations (because of bad digital beam forming, 
network connection issues, or other reasons). We therefore incor¬ 
porated conservative filtering steps into the pipeline, to identify 
and flag stations performing well outside of the normal bounds. 
The most important step considers the statistics of each station, 
and flags those that have an exceptionally large number of base¬ 
lines with high measured noise. Before LOFAR’s digital beams 
were well controlled, this step primarily removed stations with 
poorly focused beam responses. 


3.4. Target Imaging Pipeline 

The TIP stage of the MSSS pipeline combines the flagged, 
demixed, and flux-calibrated target snapshot observations to 
generate an initial image of the field. The TIP also optionally 
performs a self-calibration major cycle. Because it is decoupled 
from the pre-processing part of the pipeline, it can be run in an 
asynchronous manner with respect to the observational part of 
MSSS. 

As a first step, the individual 2 MHz band snapshots are 
combined, resulting in 8 visibility data sets per LBA or HBA 
field. These bands are treated separately throughout the TIP. The 
phase-only, direction-independent calibration is achieved using 
a VLSS-based sky model but taking the station sensitivity pat¬ 
terns into account. Example gain phases are shown in Fig.[^ The 
phases are shown for two representative stations, CS302 (about 
2 km southwest of the central group of six stations, collectively 
called the “superterp”) and RS306 (about 15 km west of the su- 
perterp). These are shown as a function of time, with one set 
of phases for each of the 8 frequency bands, and displayed here 


with an arbitrary offset for visual clarity. Phase calibration is per¬ 
formed such that an independent solution is produced for each 
2 second timestep. 


3.4.1. Imaging MSSS data 

Imaging MSSS data i s performed with the LOFAR imager, 
called the awimager ( |Tasse et~ar] 201 3|). It is based on the 
CASA imager including w-projection ( Cornwell et al.|2008] ), and 
also includes a LOFAR -specific implementation of A-projection 
( Bhatnagar et al.|2008) ) that treats the dipole and station response 
as time- and direction-dependent, full polarization terms in the 
measurement equation (e.g.^ Hamaker et al.|1996] ) during imag¬ 


ing and deconvolution. See |Tasse et al.| ( |2013| ) for the results of 
simulations demonstrating the fidelity of the imaging step. 

Our implementation of the LOFAR beam includes three lay¬ 
ers. First, the LOFAR element beam is modelled through a full 
electromagnetic (EM) simulation (not including mutual cou¬ 
pling) and implemented in our software as a polynomial fit in 
elevation, azimuth, and frequency. The HBA analog tile beam is 
also included as the direct Fourier transform (DFT) of the dipole 
positions within a tile, and rotated to account for the orientation 
of each station (as described by [van Haarlem et al.]|2Q13| the 
layout of the dipoles in each HBA station is rotated to reduce 
the sensitivity to bright far off-axis sources). Finally, the digital 
station beam for both LBA and HBA is calculated using a DFT 
of the dipole (LBA) or tile (HBA) positions in each individual 
station. Missing dipoles (LBA) and tiles (HBA) are indicated as 
such in the visibility data sets recorded by the correlator, and 
left out of the beam prediction during calculation. The digital 
station beams have been observationally mapped using the pro¬ 
cedure described by van Haarlem et al.| ( |2013 ) and found to be 
in qualitative agreement with the predictions of the beam model. 
For MSSS data, the beam is applied in the awimager such that 
it is considered to be constant within frequency blocks of width 
100 kHz, and time blocks of 5 min. 

Our initial imaging run per field incorporates projected base¬ 
lines shorter than 2kA. For fields at declination S < 35°, we 
leave out baselines shorter than 100 T, which we found empir¬ 
ically to provide a smoother background. This imaging run is 
performed with a simple, shallow deconvolution strategy (using 
2500 CLEAN iterations). After rec alibrat ing the data on the basis 
of this first imaging round (see § |3.4.3| ), we update the imaging 
parameters in preparation for final catalog creation: the maxi¬ 
mum baseline length is increased to 3 kA. After an initial shal¬ 
low deconvolution, we create a mask on the basis o f a so urce 
detection round performed in the way described in § 3.4. 2| and 
subsequently perform a deep deconvolution using the mask to 
limit the locations of CLEAN components. The final deconvolu¬ 
tion is limited by reaching a cutoff of 0.5 cr instead of by lim¬ 
iting the n umber of com ponents. All imaging steps use Briggs 
weighting ( Briggs|1995| ) with a robust parameter of 0. Final im- 
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Fig. 8. Gain phases determined by the TIP (see § |3.4| ) for observations of the HBA and LB A fields H229+70 and L227+69 (see 
Table[^. Gain phases are displayed for HBA (top row of panels) and LBA (bottom two panels), referenced to CS002HBA0 and 
CS002LBA, respectively. (CS002 is one of the six superterp stations.) The phases are shifted vertically for display purposes. They 
correspond to the (0,0) element of the station gain matrix for CS302 (top left and middle) and RS306 (top right and bottom). The 
gaps between snapshots are compressed for display purposes, and correspond to 4 hours for HBA and 45 minutes for LBA. 


age products are produced by mosaicing individual pointings in 
each band using the standard inverse variance weighting tech¬ 
nique and making use of the predicted effective primary beam 
images that are produced as a standard output of the awimager. 


We note that we expect a low impact of CLEAN bias (i.e., the 
reduction in recovered flux density of real sources due to de- 
con volution of sidelobes inadvertently ide ntified as true sources; 
e.g. [Becker et al.|1995}|Cohen et al.|2007) ). First, the snapshot uv 
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coverage is excellent for imaging at the low angular resolution 
that we make use of in MSSS (see Figures and |^, meaning 
that the synthesized beam has low sidelobe levels (rms sidelobe 
levels of < 1% in both HBA and LBA, with isolated maxima 
of 12 - 13% in the HBA and isolated maxima very close to the 
main lobe of 16 - 28% in the LBA). Second, the masked decon¬ 
volution that we employ ensures that we only apply CLEAN to 
real sources. Since sidelobes shift position at different frequen¬ 
cies, our multi-frequency source detection mitigates the impact 
of sidelobes in any individual band. We characterize the CLEAN 


bias present in MSSS-HBA data in § 5.2 


A unique benefit of incorporating a time-, frequency-, and 
direction-dependent term in the imaging and deconvolution step 
through the A-projection algorithm is the ability to directly in¬ 
corporate ionospheric corrections. We imp lement this procedure 
for our LBA data, as described in § 3.4. 3[ For the application of 
the ionospheric correction we increase the time resolution to 10 s 
in order to capture the rapidly variable phases. 


3.4.2. Source finding 

The source finding is performed using two complementary 
source extraction software packages: PyBDS]V0and PyS^ 

Both packages allow calculation of rms and mean images 
and identification of sources in radio maps through the use either 
of a False Detection Rate (FDR) method ( [Hopkins et al.|2002| ), 
or of a threshold technique that locates islands of emission above 
some multiple of the noise in the image. Gaussians are then fit¬ 
ted to each island for accurate measurement of source properties. 
Both PyBDSM and PySE allow the simultaneous use of sepa¬ 
rate detection and analysis image files for island definition and 
Gaussian fitting respectively. We make use of this functionality 
to produce a reliable catalog that retains the full multifrequency 
information provided by the data. 

While PySE has been designed as a source finding tool in¬ 
tended primarily for use by the LOEAR Transients Key Project 
and is therefore conceived for the detection of unresolved 
sources, PyBDSM has been programmed for the more general 
search of both compact and diffuse sources. PyBDSM thus al¬ 
lows fitting of multiple Gaussians to each island, or grouping 
of nearby Gaussians within an island into sources. In addition 
a PyBDSM module is available to decompose the residual im¬ 
age resulting from the normal fitting of Gaussians into wavelet 
images of various scales. This step is useful for automatic de¬ 
tection of diffuse sources. In the following, we therefore com¬ 
bine PyBDSM and PySE results only in the case of point-like 
sources. Our source finding strategy consists of running both 
PyBDSM and P ySE separately and then combining the results 
as described in § |4.2[ Eight different joint PyBDSM and PySE 
catalogs are therefore initially produced for each of the 8-band 
maps in both HBA and LBA, and these are subsequently merged 
to give the final multi-frequency catalog. 

We define two thresholds for the islands: one to determine 
the region within which source fitting is done, and another such 
that only islands with peaks above the threshold are used. Eor 
MSSS images, we set these two thresholds to 5cr and Ter, re¬ 
spectively. In addition, similarly to what is extensively done in 
the visible domain (see e.g. jSzalay et al.||1999| ), we use an 8- 
band combined image (see below for a description) and each 
single band image as detection and analysis images respectively. 
The use of a combined image for island definition optimizes sen¬ 


sitivity to faint sources. Since the significant islands are iden¬ 
tified using a single image, this procedure also alleviates the 
task of matching the eight single-band catalogs. Note however 
that the Gaussian fitting is performed on each single-band im¬ 
age independently. This results in possible differences between 
the central position of each source as a function of frequency. 
Therefore the re sultin g multi-band catalogs need to be matched 
as described in § |4.2| 

The 8-band mosaic images on which we run the source finder 
tools are produced using the same i/v-range and convolved to 
a common resolution using the miriad ( Sault et ar]|1995 ) task 
CONVOL. We do this both for primary-beam (PB) and non- 
primary-beam (NPB) corrected images. The latter are used to 
produce the combined mosaic image, obtained by performing an 
inverse-variance weighted average of the 8-band NPB corrected 
mosaic images (i.e., the weight per image is w/ = 1 ja] where ct/ 
is the rms of image /). This combined image is used as the detec¬ 
tion image, while the individual 8-band PB corrected images are 
used as analysis images. In this way, we avoid fake detections in 
the image borders (which may be caused by the increase of the 
rms related to PB correction), while obtaining properly corrected 
flux densities per source in the output of the source finders. 

Eor the combined catalog description and additional detailed 
information regarding the method by which it is produced, see 


3.4.3. Calibration stability and ionospheric corrections 

To produce the final survey output, we run the TIP twice. The 
first time incorporates the VLSS-based sky model as described 
in § |3.4| Eor the second pass, we repeat the phase calibration 
using a sky model created from the first-pass calibration and 
imaging round. On the basis of this new phase calibration, we 
generate a new set of images (as described in § |3.4.1| ) and source 
catalogs. The intention of this step is to minimize the effect of 
any spurious sources that may be present in the initial sky model 
and to ensure that the MSSS catalog is based on an internally 
consistent calibration cycle. This procedure has been followed 
for the representative data set considered in §[^ 

Eor the low frequency and large field of view intrinsic to 
LOEAR observations in the LBA band, ionospheric effects are 
strong even when imaging at the modest 2' resolution uti¬ 
lized for MSSS. The clearest effect in the image plane is the 
presence of spiky artefacts around the brightest sources. To deal 
with this we have implemented a scheme si milar to the Source 
Peeling and Atmospheric Modeling (SPAM; [rntema et al.|2009| ) 
approach that has been successfully used for VLA and GMRT 
data. We briefly summarize the procedure here; a full description 
will be provided in a forthcoming publication. The sky model is 
divided into approximately 30 source groups, and each group 
is used to derive a phase solution at each frequency. These 30 
source groups cover the entire field of view visible at 30 MHz. 
To trace the frequency behavior of the calibration phases in those 
directions at the highest frequencies, where the field of view is 
much smaller, we utilize the simultaneously observed flanking 
beams. Thus ionospheric phases are only available across the 
full field of view for the central field in a multiplexed observa¬ 
tion like MSSS. Eor the dataset considered in §[^ only the central 
field (L227-h69) can be processed in this manner. 

The frequency dependent phases include two terms: a non- 
dispersive clock delay ternj^ and a dispersive ionospheric delay 


^ http: //tinyurl. com/PyBDSM-doc ^ lqFAR stations only share a single clock within the core area; the 

http: //docs. transientskp. org/tkp/master/tools/pyse. htnrimote stations are on independent clocks. 
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term. The clock term is constant in all directions at a given time, 
while the dispersive term is direction-dependent. To isolate the 
ionospheric term we subtract the phases determined in one di¬ 
rection from those in all other directions, leaving a differential 
ionospheric phase term per direction and per station. By con¬ 
sidering the phase of each station in each direction as a single 
“pierce point” through a single thin-layer ionosphere, these val¬ 
ues are used to fit a “screen” of ionospheric total electron content 
(TEC) at a particular height. The difference between data and 
fitted screen is used to identify an optimized TEC screen height, 
typically around 100-200 km. The TEC screen can be used as an 
input to the awimager to correct the phase distortions across the 
field of view, as a function of both frequency and time. 


4. Standard data products 

The primary output of MSSS is a catalog containing positions 
and Stokes I flux densities for all confidently detected sources, 
as well as extents and orientations for resolved sources. Spectral 
behavior of the sources is also provided. The MSSS catalog is 
stored in the LOEAR GSM database, implemented in a fully 
VO-compliant system based on the Data Center Helper Suite 
(DaCHS; |Demleitner|20T^ . 

Processed visibilities will be stored in the LOEAR long term 
archive (ETA) for future reprocessing, for example by science 
groups that wish to reprocess the data to search for polarization, 
or to perform long-baseline imaging (see §|^. In addition, raw 
visibilities are always stored in the ETA immediately after obser¬ 
vation. At the conclusion of the TIP, the images are imported to a 
postage-stamp serveij^ which allows inspection of images with 
catalog overlays and multi-frequency source spectrum pop-up 
plots, as well as providing direct-download links to PITS files. 

Pollowing verification of the survey results, the survey out¬ 
put (images and catalogs) will become fully public. 


4.1. MSSS Images 

The MSSS image products will be released in FITS format, and 
will c onsist of mosaics corresponding to the fields specified in 
§ 2.4 Each mosaic consists of sixteen 2 MHz bandwidth images 
at the central frequencies listed in Table and two 16 MHz full 
bandwidth images, one for LBA and the other for HBA. 


4.2. MSSS Catalog 

The MSSS catalog is generated in several steps. Pirst, the PySE 
and PyBDSM source finders are run on all individual HBA and 
LBA frequency b ands u sing the combined maps as detection im¬ 
ages (see Section [3A2 for details). When deriving source posi¬ 
tions, sizes and fiuxes, both PyBDSM and PySE take into ac¬ 
count fitting errors caus ed by correlated map noise following the 
approach suggested by Condon ( |1997| ). Ionospheric phase cal¬ 
ibration is the other largest contributor to the positional uncer¬ 
tainty of fitted sources ( [Cohen et al.||2007| ). We describe in the 
following how these errors have been taken into account when 
producing the multi-band final catalog. 

We firstly work on the HBA and LBA catalogs separately. 
Por each source finder, the detected sources are associated across 
the eight individual frequency bands in order to generate a con¬ 
catenated multi-frequency source catalog. The association is per¬ 
formed by calculating the angular distance between sources in 


^ MSSS data are being hosted at http: //vo. astron. nl where the 
subset presented in §[^has already been made available. 


any pair of individual frequency band catalogs. A match is de¬ 
termined to be positive if each element of the pair is mutually 
the nearest to its counterpart from the other catalog. This cri¬ 
terion ensures that sources have exclusive pairing, as opposed 
to possible multiple associations. An additional threshold is ap¬ 
plied in order to reject sources that are too distant and might be 

spuriously associated; hence we require distance < 3 ^Jcr^ cr^, 

where cr 12 are the fitted positional uncertainties in the first and 
second catalog, respectively. Note that, in this step, we do not 
take into account calibration related errors, since in each of 
the two segments (HBA and LBA) the 8 bands have been ob¬ 
served simultaneously. We chose a threshold of 3cr after verify¬ 
ing that the curve of growth of positive matches started plateau- 
ing around this level. After the source association is completed, 
we calculate the position of each source - RA and Dec separately 
- as the weighted average position among the frequency bands 
in which it was detected taking into account the positional uncer¬ 
tainties. The uncertainties associated with the average positions 
are calculated by propagating the errors accordingly. 

Eollowing the source association that has so far been per¬ 
formed only within the PySE and the PyBDSM catalogs, sources 
detected with each source finder are also cross-matched in or¬ 
der to generate the final catalog. The same pairing process as 
explained above is used for this step. However, in order to pre¬ 
serve consistency in the final catalog, we did not attempt to cal¬ 
culate average values for the various reported fields. Rather, we 
chose the PyBDSM fields to prevail over those from PySE when 
both values existed. Hence, the reported IDs and positions are 
taken from PyBDSM unless a source was only detected by PySE. 
Similarly, for a given frequency band, the reported fiux den¬ 
sity properties are those found by PyBDSM unless the source 
was only detected in this particular frequency band by PySE. 
Each frequency band possesses a field SFFLAG^^^ that indicates 
which, if not none or both, source finder the source was detected 
with at frequency nnn MHz. By considering the results from both 
source finders together in this way, we can add confidence to 
the reality of individual detections; thus we recommend using 
sources detected by only one source finder with caution. 


In order to produce the final source list and associate the 
HBA and LBA catalogs, we need to incorporate an estimate 
of the eff ect on source positi ons caused by calibration errors. 
Eollowing [Cohen et al. ( 2007[ ), we compare our HBA and LBA 
catalo gs to the NRAO/VLA Sky Survey (NVSS; [Condon et 
1998[ ) catalog, which has higher resolution and S/N, as well as 


significantly lower calibration errors due to the higher observing 
frequency (1.4 GHz). At each of the two frequencies we iden¬ 
tify sources with i) a detection level of at least 30cr, ii) a sin¬ 
gle brighj^NVSS counterpart within 1.5 the beamsize, and iii) 
a fitted major axis less than 1.5 the beamsize. Eor example, in 
the case of the LBA data presented in §|^ we derive an average 
offset of ARAmean = O'.'18 and ADeCmean = O'.'03 with associ¬ 
ated rms deviations from the mean values of ARAnns = 1''59 
and ADeCrms = O'.'24. Eor the HBA data in §[^ where no di¬ 
rection dependent calibration was applied, the coordinates dis¬ 
play larger offsets: ARAmean = 2'.'18 and ADeCmean = -O'.'79 
with rms deviations from the mean values of ARAnns = 2'.'92 
and ADeCrms = 2'.'45. During the production of the final multi¬ 
band (16-frequency) catalog, these calibration errors are added 
in quadrature to the positional uncertainties of the fitted sources. 
The HBA and LBA association is subsequently performed as de- 


^ Peak flux higher than 50 mJy beam ^. 
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scribed before, but taking into account both fitting and calibra¬ 
tion position errors. 

Finally, we perform a post-concatenation analysis in order to 
determine the spectral properties for each source. For this step, 
only the PyBDSM fluxes are used - again this is done in order to 
ensure consistency since the flux scale between the two source 
finders may suffer from biases. The spectrum of each source is 
fitted with the functional form (see also |Scaife & Heald|2012| ): 

= (6) 


Given the large number of sources to be fitted and the fact 
that the posterior distribution of the spectral fit should be well- 
behaved - this is a linear least-squares fit to a polynomial in 
log^y space - we use a Levenberg-Marquardt minimisation 
algorithm to determine the best-fit parameters and errors. 

We use a locally determined effective beamsize to de¬ 
convolve the source sizes reported by PyBDSM, and classify 
sources as extended if the deconvolved size is nonzero. Since 
most sources are unresolved at this resolution, this procedure al¬ 
lows us to mitigate the effect of ionospheric smearing. 

The 214 MSSS mosaics overlap at their edges to ensure that 
all sources are reliably imaged and cataloged. This means that 
many sources are identified in more than one mosaic. After creat¬ 
ing a catalog from each mosaic as described above, we filter the 
multiply cataloged sources to remove duplicate entries. First, we 
look for sources that have the same source ID, which identifies 
those having the same coordinates at arcsecond precision in both 
RA and Dec. Since the mosaics are formed from the same im¬ 
ages, most sources that are present in more than one mosaic are 
identified at the same coordinates. However, small differences 
sometimes arise because the local backgrounds and noise lev¬ 
els calculated by the source finders differ slightly in neighboring 
mosaics. Therefore, an extra sifting is performed by identify¬ 
ing the nearest neighbour of each cataloged source. Such neigh¬ 
bors are matched and removed if they were found in different 
mosaics, and if their separation is less that 45" (substantially 
smaller than the resolution, but large enough to identify match¬ 
ing source pairs even when their central position is less precise). 

The MSSS catalog has the following columns. First, a set of 
parameters that are common for both the point source catalog 
and extended source catalog: 

ID Source ID, formed as ‘TMSSS ^hhmmss + 

ddmmss'' using the lAU convention. In the cat¬ 
alog presented in §|^ we instead use the source 
ID prefix “MSSSVF” to distinguish it from the 
forthcoming full MSSS catalog. 

RA Source J2000 Right Ascension, in decimal de¬ 

grees. 

DEC Source J2000 Declination, in decimal degrees. 

e_RA Error in source J2000 Right Ascension, in sec¬ 

onds of time. This is a formal error based on 
the source position fit. 

e_DEC Error in source J2000 Declination, in arcsec- 

onds. This is a formal error based on the source 
position fit. 

e_RA_sys Eull error in source J2000 Right Ascension, in 

seconds of time. This includes both the formal 
error based on the source position fit and a sys¬ 
tematic positional error term. 

e_DEC_sys Eull error in source J2000 Declination, in arc- 
seconds. This includes both the formal error 
based on the source position fit and a system¬ 
atic positional error term. 


SYYLkQnnn 


Sirvtnnn 

eJS>i.rvtnnn 

e_Spk/t^^ 

A0_LBA 

e_A(S)_LBA 

A1_LBA 

e_Al_LBA 

A(S)_HBA 

e_A(S)JlBA 

A1_HBA 

e_AlJlBA 

AO 

e_A0 

Al 

e_Al 

NDET 

NDET_LBA 

NDET_HBA 

NUNRES 

NUNRES_LBA 

NUNRES_HBA 

MOS_ID 

CAL_ID_LBA 

CAL_IDJ1BA 


Elag indicating which source finder identified 
the source at nnn MHz (0 means it was de¬ 
tected in both; 1 means it was detected only in 
PyBDSM; 2 means only in PySE; 3 means no 
detection). Note that if the source was identi¬ 
fied by both source finders, the reported flux 
density values are those of PyBDSM. 

Source integrated flux density at nnn MHz, in 
Jy. 

Error in source integrated flux density at nnn 
MHz, in Jy. 

Source peak flux density at nnn MHz, in 
Jy beam“^. 

Error in source peak flux density at nnn MHz, 
in Jybeam"^ 

Spectral flux density at 150 MHz, Aq in 
Equation!^ Derived from LB A values only. 
Error in spectral flux density at 150 MHz. 
Derived from LBA values only. 

Spectral index, Ai in Equation]^ Derived from 
LBA values only. 

Error in spectral index. Derived from LBA val¬ 
ues only. 

Spectral flux density at 150 MHz, Aq in 
Equation!^ Derived from HBA values only. 
Error in spectral flux density at 150 MHz. 
Derived from HBA values onW 
Spectral index, Ai in Equationl^ Derived from 
HBA values only. 

Error in spectral index. Derived from HBA 
values only. 

Spectral flux density at 150 MHz, Aq in 
Equation!^ 

Error in spectral flux density at 150 MHz. 
Spectral index, Ai in Equation 
Error in spectral index. 

Number of bands in which the source was de¬ 
tected. 

Number of LBA bands in which the source 
was detected. 

Number of HBA bands in which the source 
was detected. 

Number of bands in which the source is unre¬ 
solved. 

Number of LBA bands in which the source is 
unresolved. 

Number of HBA bands in which the source is 
unresolved. 

Mosaic name from which the source was ex¬ 
tracted. 

Primary flux calibrator(s) used to set the LBA 
flux density scale in the vicinity of the source. 
Primary flux calibrator(s) used to set the HBA 
flux density scale in the vicinity of the source. 


A set of parameters that are only included for the extended 
source catalog: 


MAJAX^^/t 

MINAX^^^ 

Vknnn 

e_MAJAX^/t/t 


Major axis of fitted ellipse at nnn MHz, in dec¬ 
imal degrees. 

Minor axis of fitted ellipse at nnn MHz, in dec¬ 
imal degrees. 

Position angle of fitted ellipse at nnn MHz, in 
decimal degrees. 

Error in major axis, in decimal degrees. 
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e_MINAX/t^^ Error in minor axis, in decimal degrees. 
e_PA^^^ Error in position angle, in decimal degrees. 

This catalog has 108 columns for the point sources, and 108+ 
96 = 204 columns for extended sources. Assuming that there 
will be 10^ sources (see §i with about 5% of those extended, 
then that will lead to about 11.3 million cataloged data values. 


5. Initial MSSS imaging results 


MSSS observations are typically performed for several hours per 
week within the rest of the LOEAR schedule. The HBA seg¬ 
ment has been fully observed and initial images have been cre¬ 
ated. While the survey calibration processing is still ongoing, 
we highlight a representative portion of the sky to illustrate the 
output that will be forthcoming for the full northern sky. Eor 
this we selected a 100 square degree patch of sky that has sub¬ 
sequently been used repeatedly to test imaging pipeline perfor¬ 
mance. The field was randomly selected and we refer to it here 
as the MSSS Verification Eield (MVE). It includes a few mod¬ 
erately bright point-like sources but no bright, complicated 3C 
or 4C sources, and is otherwise distant from the troublesome 
A-team sources. The Galactic contribution to the sky brightness 
is relatively unimportant in this direction (Galactic coordinates 
(/ = 108°, b = 44°) at the center of the MVE). 

The mosaics in this region are created from 9 LBA fields and 
32 HBA fields, as illustrated in Eig.[^ These fields were not all 
observed at the same time or even on the same day. The observ¬ 
ing summary is listed in Table The primary fiux calibrator for 
fields H244+70 and H245+75 was 3C380, and the primary fiux 
calibrator for all other fields was 3C295. After convolution to a 
common beam size, the effective resolution is 108" in HBA, and 
166" in LBA. Images of the frequency-averaged LBA field and 
HBA mosaic are shown in Eigs.p^andpJ^ respectively. 

Due to the ionospheric processing scheme described in 
§ |3.4.3| only the central LBA field (L227+69) is used to gen¬ 
erate the LBA part of the MVE catalog. Moreover, only 7 of the 
9 LBA snapshots were used due to particularly bad ionospheric 
quality during the first 2 snapshots (see Eig. [^. The HBA por¬ 
tion of the catalog was produced based on the combination of all 
fields listed in TableAll told, 1234 unique sources were iden¬ 
tified within the 100 square degrees of the MVE (307 in LBA, 
1234 in HBA). A simple projection to the full MSSS survey area 
would suggest that approximately 250 000 sources will be found, 
but taking into account reduced sensitivity at low declination and 
Galactic latitude, as well as poor image quality near extremely 
bright sources, we expect to recover between 150000-200000 
sources in the full MSSS catalog. 

We present the spectral index distribution determined on the 
basis of the MVE catalog in Eig. [T^ We consider all sources 
with a cataloged AO value greater than 200 mJy (628 out of 
1209 sources). Considering the spectral index determined us¬ 
ing all cataloged frequencies (recorded in column Al), the mean 
and median values are -0.60 and -0.66, respectively. These are 
somewhat more shallow than spectral indices determined for 
the same sources from HBA fluxes alone (mean and median of 
-0.70 and -0.77, respectively). Part of this may be due to spec¬ 
tral turnovers at LBA frequencies or cosmic ray energy loss pro¬ 
cesses, although some part is likely due to measurement errors. 
We note that the HBA-only spectral indices have a systematic 
error due to beam effects not incorporated in the beam model de¬ 
scribed in § |3.4.1| The error is estimated at the level of 0.05+0.22 
for the MVE region on the basis of more recent electromagnetic 
simulations of the HBA stations including the effects of mutual 




HBA spectral index (A1_HBA) 


Fig. 12. Spectral index histograms for the MVE catalog sources 
with A0_HBA> 0.2 Jy, based on the cataloged Al values (top), 
and the same but for HBA values only, catalogued as AIJHBA 
{bottom). 


coupling. Because these simulations are still under development, 
and the correction at the high elevation of the MVE snapshots is 
expected to be very small, we have not adjusted the HBA spectral 
indices presented in the bottom panel of Eig. We will fully 
address this issue during the creation of the full MSSS catalog. 

5.1. Completeness and false detection rate 

We determined the completeness of the MVE portion of the sur¬ 
vey in the standard way through the use of Monte Carlo sim¬ 
ulations, injecting simulated point sources into residual images 
from the survey and attempting to recover them with the source 
finder. Note that this approach only considers systematic issues 
related to source identification and characterization in the im¬ 
age plane. We used PyBDSM for this process; results with PySE 
would be expected to be very similar. The only complication 
arises from the fact that averaged images (from all LBA and 
HBA bands) are used as detection images in the cataloging pro¬ 
cess. Therefore, as a major step in finding the completeness in 
individual bands, we replicated the full procedure used to gener¬ 
ate and apply detection images during the source finding process. 


13 














































Heald et al.: LOFAR MSSS 




RA (J2000) RA (J2000) 

Fig. 9. Mosaic layout of the MSSS Verification Field (MVF). Left: The nine LB A fields making up the MVF, overlaid on the 10° x 10° 
mosaic field (central white square). The gray border around the mosaic area is a guard area used to ensure that sufficient edge fields 
are included and to ensure fiat sensitivity within the mosaic. Right: The 32 HBA fields making up the MVF, overlaid on the same 
central mosaic field. The guard border is smaller in proportion to the HBA field diameter. 
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Fig. 10. LBA full-bandwidth central field of the MVF, after ionospheric correction and displayed here without primary beam cor¬ 
rection for clarity. The noise level is 39 mJy beam“\ and the synthesized beam is 166". The colorbar units are Jy beam“^ Diamonds 
mark the positions of cataloged VLSSr sources. 
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Fig. 11. HBA full-bandwidth mosaic of the MVF, displayed here without primary beam correction for clarity. The noise level is 
5mJybeam“\ and the synthesized beam is 108". The colorbar units are Jybeam“^ Diamonds mark the positions of cataloged 
VLSSr sources. 


Table 3. MSSS Verification Field Observing Log 


Field 

Date 

Field 

Date 

Field 

Date 

Field 

Date 

Field 

Date 

L212+63 

2013 Apr 8 

H214+63 

2013 Apr 21 

H231+65 

2013 Apr 21 

H215+70 

2013 Apr 21 

H228+73 

2013 Feb 15 

L225+63 

2013 Mar 18 

H220+63 

2013 Apr 21 

H237+65 

2013 Apr 21 

H222+70 

2013 Feb 15 

H236+73 

2013 Feb 15 

L238+63 

2013 Apr 8 

H225+63 

2013 Apr 21 

H212+68 

2013 Apr 21 

H229+70 

2013 Feb 15 

H244+73 

2013 Apr 21 

L211+69 

2013 Mar 18 

H230+63 

2013 Apr 21 

H218+68 

2013 Feb 15 

H236+70 

2013 Feb 15 

H208+75 

2013 Apr 21 

L227+69 

2013 Mar 18 

H236+63 

2013 Apr 21 

H224+68 

2013 Feb 15 

H244+70 

2013 Feb 22 

H217+75 

2013 Apr 21 

L243+69 

2013 Mar 18 

H214+65 

2013 Apr 21 

H231+68 

2013 Feb 15 

H204+73 

2013 Apr 21 

H226+75 

2013 Apr 21 

L201+75 

2013 Apr 8 

H220+65 

2013 Feb 15 

H237+68 

2013 Feb 15 

H212+73 

2013 Apr 21 

H235+75 

2013 Apr 21 

L222+75 

L244+75 

2013 Mar 18 
2013 Apr 8 

H226+65 

2013 Feb 15 

H208+70 

2013 Apr 21 

H220+73 

2013 Feb 15 

H245+75 

2013 Feb 22 


In detail, we generated residual maps from the actual images 
used for the cataloging, after first removing any detected sources 
with PyBDSM: this ensures that the noise and its spatial distri¬ 
bution are consistent with that of the real data. We drew sim¬ 
ulated sources from a power-law flux density distribution with 
with fixed upper and lower fiux densities chosen to 
span the range of observed fiux densities in the survey. For the 
multi-band analysis we additionally drew source spectral indices 
from a Gaussian distribution with mean -0.7 and standard devi¬ 
ation 0.35, and considered the fiux density to be at a reference 
frequency in the middle of the band (135 MHz for HBA, 50 MHz 
for LB A). A suitable number of simulated sources were then 
added at random positions to the residual maps for each band. 


For the multi-band analysis, we constructed a detection image 
by averaging the individual bands (in the case of the HBA data, 
this was done by using residual images without beam correction, 
taking account of the beam correction factor by scaling the input 
fiuxes) in order to mimic as closely as possible the process used 
in cataloging. Finally, PyBDSM was used to attempt to recover 
the simulated sources from the resulting image, taking care to 
use exactly the same parameters as applied in the cataloging: a 
source was deemed to have been recovered if PyBDSM detected 
a source within 1 arcmin of the input position with a fiux density 
that matched to within lOcr, where cr is the fiux density error re¬ 
turned by PyBDSM. (In practice, positions normally matched to 
within less than one pixel.) This process gave a list of matched 
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and unmatched sources, which, after repeating several times to 
improve the statistics, could be used to estimate the survey com¬ 
pleteness. 

The results of the completeness simulations for LB A and 
HBA are shown in Fig. 13 This figure shows cumulative com¬ 
pleteness curves, i.e. it gives the completeness for sources above 
the indicated fiux density limit. We see that the HBA catalog is 
expected to be 90% complete above 100 mJy, and 99% com¬ 
plete by around 200 mJy at a mid-band reference frequency 
of 135 MHz; the cataloging process using the averaged detec¬ 
tion image gives, as expected, a catalog that is roughly a factor 
Vs more sensitive than those derived for the individual bands. 
The improvement in sensitivity realized after frequency averag¬ 
ing suggests that the MSSS-HBA images have not yet reached 
the classical confusion limit (cf. § |2.1| ). The LB A catalog has a 
much higher completeness threshold of 0.55 Jy (90% complete) 
or 0.80 Jy (99% complete) at a reference frequency of 50 MHz. 

We emphasise that these completeness curves are for point 
sources only (though at the resolution used that includes nearly 
all real sources), that the process assumes that the beam is well 
modelled as a Gaussian, and that residual images are free from 
real structure, which are good assumptions for the HBA images 
but much less so for the LBA. We also assume that there are 
no relative fiux scale offsets within the LBA and HBA bands. 
Any departures from these assumptions will tend to make the 
real catalogs less complete than indicated by the completeness 
curves. At present we regard the HBA curves as reliable, but the 
LBA curves should be taken as indicative only. We note that the 
simulation used for our completeness estimates only considers 
effects present in the image plane. During the development of 
the awimager, simulations of ideal point sources were used to 
demonstrate e xcellent image pla ne recovery of objects added to 
the visibilities ( |Tasse et al.|2013| ), so we do not expect a substan¬ 
tial effect on the completeness due to issues in the imaging soft¬ 
ware that we use. Imperfections in our calibration solutions and 
beam model may negatively impact completeness, but a more 
detailed simulation to address those effects is beyond the scope 
of this paper. A forthcoming paper will present the MSSS cata¬ 
log over the full survey area, and with that much larger statistical 
sample some of these issues may be more effectively addressed. 

Finally, a by-product of this simulation process is a test 
of the reliability of source fiux densities as extracted with 
PyBDSM. We show a representative plot in Fig. [T^ It can be 
seen that PyBDSM recovers the fiux density very accurately. A 
few sources at low fiux densities are found to have significantly 
(a factor of a few) high fiux densities relative to the input values: 
we attribute this to confusion (i.e. there is some overlap with a 
nearby bright source which is not completely deconvolved by 
PyBDSM). As the right-hand panel of Fig. shows, however, 
such sources are a very small fraction of the total. 

Having addressed the completeness of the MVF catalog, 
we now proceed to assess the possibility of falsely detected 
sources being included in the catalog. To that end we have cross- 
matched the MVF catalog with existing radio surveys, primar¬ 
ily the deeper 1400 MHz NVSS catalog. Of the 1209 sources in 
the MVF catalog, we find that all but 8 are associated with an 
NVSS source. That result is based on searching for counterparts 
within the MSSS resolution element and a visual comparison to 
recover associations with components of extended sources. The 
8 remaining MSSS sources were then compared with VLSSr and 
the 151-MHz 7C survey ( [Hales et ST]|2007| see § |5.3.1| ). Three 
correspond to VLSSr sources, but none were found in 7C. Thus 
there are 5 MVF sources that have no associated previously cata¬ 
loged radio source. All have fiux densities between 50 - 200 mJy 


and could be steep spectrum sources below the detection thresh¬ 
olds of both NVSS and VLSSr (a spectral index of a ^ -2 would 
be sufficient), but we nevertheless consider them to represent an 
upper limit to the MSSS false detection rate (FDR), < 0.4%. 


5.2. CLEAN bias 


Here we be gin to characterize the CLEAN bias (see, e. g., [Becker 
et al. 1995\ present in the MSSS catalog. As noted in § 3.4.1[ we 
anticipate a small bias because of the excellent uv coverage and 
our careful deconvolution procedure. To characterize the effect 
in MSSS data, we have performed a simulation adding artificial 
sources into the visibilities, followed by an imaging sequence 
carried out in the same way as for the actual data. We started 
with the visibility data from all 8 HBA bands of a representa¬ 
tive field in the MVF region (H229+70). Real sources that are 
characterized in the MVF catalog were subtracted in the vis¬ 
ibility plane. Artificial point sources were then drawn from a 
p opul ation similar to the one used for the completeness study in 


5.1 


W here spanning fiux densities between 30mJy 

and 30 Jy). These were added into the visibilities, incorporating 
the LOFAR beam model. Next, the visibility data were imaged 
as described in §[3.4. 1[ In particular, an initial shallow clean of 


each band was used to generate a stacked full-bandwidth im¬ 
age that defines a mask for a subsequent deep deconvolution of 
each separate band. We used the awimager for the imaging steps 
including the full LOFAR beam model. We also verified our re¬ 
sults by injecting the same sources into the visibilities without 
incorporating the LOFAR beam model and imaging the result¬ 
ing visibilities in casa with the same settings as were used in 
awimager (e.g. the same CLEAN mask). 


Pixel values were drawn from the known locations of in¬ 
jected sources and compared with the input fiux density values. 
We followed this procedure to eliminate complications due to 


but would have affected our CLEAN bias estimate had we re¬ 
lied on source finding to determine the reconstructed fiux of 
each artificial source. We find a typical CLEAN bias of lOmJy 
in each band, well below the typical noise level in each band 
30mJybeam“^) and far below the completeness threshold. 
Multiple realizations of the same artificial source population 
were utilized to assess the robustness of the estimate. We found 
that typical values cluster around lOmJy, but with substantial 
variation around the 50% level, depending on the particular re¬ 
alization of the artificial source population. For the MVF, we do 
not make a CLEAN bias correction in the catalog since its exact 
value is uncertain but very small (< 5%) above the completeness 
limit per band. In the full MSSS catalog, we will characterize 
this effect in more detail along with the completeness considera¬ 
tions described in § [5.1[ We have not characterized this effect in 
the LBA portion of the MVF catalog because it is not yet pos¬ 
sible within our direction dependent calibration procedure. It is 
expected to be a small effect for the same reasons as in the HBA 
catalog. 


completeness effects, which is characterized separately in § 5.1 


As expected, the CLEAN bias measured in MSSS-HBA is 
far lower than in other recent radio surveys. For example, in 
the original VLSS survey catalog, the bias was characterized as 
L39(T (variable with the local rms noise level), and substantially 
improved in VLSSr to 0.66cr. In the same terms, the CLEAN bias 
is 0.67O- in NVSS and L67o- in FIRST. For MSSS-HBA, we find 
a CLEAN bias of approximately 0.3cr using this simulation. 
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Fig. 13. Results of completeness simulations for LB A (left) and HBA (right). Curves show the cumulative completeness for the 
individual frequencies in the LB A and HBA and for their combination at a nominal reference frequency of 50 and 135 MHz, 
respectively. 
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Fig. 15. Comparison of 7C and MSSS cataloged sources in the 
MVF. 


5.3. Comparison with other catalogs 
5.3.1. HBA 

The HBA MSSS catalog can be com pared directly to the 151- 
MHz 7C survey ( [Hales et~ar]|2007| ). 7C has a resolution of 
70 X 70cosec(^) arcseconds at a frequency of 151 MHz, and 
covers 1.7 sr of the northern sky, including the area of the MSSS 
described here. 

We filtered the full 1C catalog to produce a catalog within 
the 100 square degrees of the field described in this paper, con- 


^ Available electronically from 

http: //WWW .mrao.cam.ac.uk/facilities/surveys/7c/ 


taining 607 1C sources. The completeness limit for 7C in this 
part of the sky is about 400 mJy, and about 248 1C sources lie 
above this completeness limit. We would expect essentially all 


of these to be detected by MSSS since, as noted in Section 5.1 


the MSSS is 99% complete to 200 mJy at 135 MHz. The rms 
noise in the 1C images is about 30 mJy, while the noise in the 
single-band HBA images is generally less than 20 mJy, and so 
in general we expect MSSS to go deeper than 1C. Indeed, after 
filtering out cataloged sources with poorly constrained positions 
and those without measured 151 MHz fiux densities, there are 
1101 MSSS sources in the area of interest, a factor of 1.8 more 
than in 7C. 


We cross-matched the 1C and MSSS catalogs by combin¬ 
ing the random (not systematic) positional errors of each pair of 
sources in quadrature, and finding the maximum-likelihood 1C 
match for each MSSS object. Bearing in mind that the distribu¬ 
tion of offsets is a Rayleigh distribution, we imposed a threshold 
on the acceptable maximum likelihood, determined by inspec¬ 
tion of the matching results, to reject spurious matches. This 
method makes optimal use of the known positional errors. The 
initial cross-match resulted in a mean positional offset of 7.6+0.8 
and -0.8 ±0.3 arcseconds in RA and Dec respectively: we re¬ 
moved this offset (which is consistent with the offset already de¬ 
termined above) and repeated the cross-matching process. A to¬ 
tal of 592 sources, or 98% of the 7C sources, matched between 
the two catalogs, and all but 4 of the 1C sources above the com¬ 
pleteness limit (98%) are matched with MSSS sources. 

It is useful to look at the positions and fiux densities of 
matched and unmatched sources (Fig. [I3- We see that most un¬ 
matched sources are faint MSSS sources, which is just a conse¬ 
quence of the fact that MSSS is more sensitive than 1C over most 
of the sky area. The very few bright unmatched 1C sources are 
almost all positionally nearby to unmatched MSSS sources of 
similar brightness, which means that these are almost certainly 
the same sources with positions that differ by too much for the 
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Input flux density (Jy) Deviation from expected flux density (a) 


Fig. 14. Comparison of input and output model source flux densities for one band of the HBA (left) and histogram of the deviation 
from the expected value (right). The green line on the left shows the line of equality, on the right the expected normal distribution if 
the fluxes are simply affected by the noise estimated from the fit. 


algorithm to consider them to match. Sources of unmodelled po¬ 
sitional uncertainty include the slightly higher resolution of the 
7C survey. 

Finally, we can consider the flux scales of the two catalogs. 
Here we consider only the subsample of the catalog above the 7C 
completeness limit of 400 mJy, since the flux densities of parts 
of the catalog where the samples are not complete will tend to be 
biased high. There is an excellent correlation (Fig.[^ between 
the flux densities from the two samples. A very few sources have 
much larger total flux densities in 7C compared to MSSS, but 
these differences are not present in the peak flux-density analy¬ 
sis, and so must arise from differences in the algorithms for fit¬ 
ting to extended sources in the two surveys. We see that there is 
a slight but clearly significant difference in both the peak and in¬ 
tegrated flux densities in the 7C and MSSS catalogs, in the sense 
that the MSSS flux densities are systematically high by about 
9% (total flux density) or 6% ( peak flux density). Both surveys 
should be tied to the flux scale of |Roger et al. ( |1973| ), so in princi¬ 
ple we would not expect this systematic difference. In practice, 
the difference is likely due to different assumptions about the 
flux density of 3C295, the reference flux calibrator for almost 
all the MSSS observations, which seems likely to have been the 
primary flux calibrator for the 7C observations as well given the 
RA of the held ( [McGilchrist et al.||1990| ). The 7C catalog uses 
the 6C flux density of this object, 89.8 Jy, whereas the 150 MHz 
flux density interpolated from the fits of |Scaife & Heald| \20\2) 
is 97.7 Jy, giving the correct direction and approximately the 
observed magnitude of the flux density discrepancy. Indeed, the 
150 MHz flux density of 3C295 appears anomalously low on 
the flux plot of Scaife & Heald, compared to, e.g., the 178 MHz 
3CRR flux density. We are therefore confident in the flux scale 
and flux recovery in the catalog: further investigation of the true 
flux density of 3C 295 at HBA frequencies would be desirable. 


5.3.2. LBA 

LBA comparisons are restricted by the relatively small num¬ 
ber of cataloged sources with low-frequency flux densities. 
However, we have compared the MSSS MVF with the 8C at 38 
MHz and the VLSS at 74 MHz in the same way as in the previ¬ 
ous section. We And that there are a similar number of MSSS and 
8C sources in the held (279 MSSS sources, 233 8C sources), but 
only around 80% of MSSS sources (above ~ 0.8 Jy at 50 MHz, 
the 99 per cent completeness level) have a counterpart in 8C, 
and similarly there are 8C sources with no MSSS counterpart. 
The 8C flux densities are higher than the MSSS ones by about 
11%, with significant scatter. By contrast, the VLSS appears to 
go slightly deeper than MSSS, but at the 99% completeness level 
nearly all (94%) of the MSSS sources have VLSS counterparts, 
with a good agreement between peak fluxes - the VLSS sources 
are brighter by about 3%, with little scatter. A detailed compari¬ 
son with these catalogs is deferred until a larger held is available. 
We note that while the MSSS data will ultimately allow a deeper 
catalog, the current images are limited by substantial ionospheric 
errors that are only partially corrected with our direction depen¬ 
dent calibration and imaging procedure. Moreover the LBA held 
is composed of just a single pointing, limiting the full sensitivity 
to the center of the 10° x 10° region considered here. We ex¬ 
pect the final LBA data products to provide substantially deeper 
levels than both the 8C and VLSSr surveys. 


6. Scientific capabiiity 

To assess how MSSS compares with other existing, ongoing, and 
future surveys, we list key parameters in TableThe survey pa¬ 
rameters illustrate that MSSS is complementary to existing cat- 
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Fig. 16. Comparison of total (left) and peak {right) flux densities for bright sources in the MSSS and 7C surveys. 


Table 4. Nominal MSSS parameters and comparison with other surveys 


Survey 

Frequency 

Sensitivity 

Resolution 

Area 

MSSS-LBA 

30-78 MHz 

< 50 mJy beam"^ 

< 150" 

20000 {6>0°) 

8C 

38 MHz 

200 - 300 mJy beam"^ 

4.5' X 4.5' csc(^) 

3 000 {6 > +60°) 

VLSS 

74 MHz 

100 mJy beam"^ 

80" 

30000 0° (^> -30°) 

MSSS-HBA 

120-170 MHz 

<10-15 mJy beam"^ 

< 120" 

20000 0° ((i> 0°) 

7C 

151 MHz 

20 mJy beam"^ 

70" X 70" csc{6) 

5 500 0° (irregular coverage) 

TGSS 

140-156 MHz 

7-9 mJy beam"^ 

20" 

32000 0° (^> -30°) 

WENSS 

330 MHz 

3.6 mJy beam"^ 

54" X 54" csc((i) 

10000 0° ((i> +30°) 

NVSS 

1400 MHz 

0.45 mJy beam"^ 

45" 

35 000 0° (6 > -40°) 


Note. Sensitivity and resolution values for the MSSS survey components are upper limits corresponding to images produced with baselines 
shorter than 3 kd. Longer baselines are included in the observations as a matter of course, enabling reprocessing toward the production of an 

updated, higher angular resolution catalog. 


alogs, with potential for even better image quality with followup 
reprocessing steps. 

The values for sensitivity and resolution are shown graphi¬ 
cally in Figure [T^ The sensitivity panel shows that sources with 
typical spectral indices will be detected in both WENSS and 
MSSS. Steep spectrum sources detected in NVSS will also be 
recovered in MSSS. Therefore, the sensitivity of MSSS make 
it interesting in terms of comparison with other all-sky surveys. 
The study of steep spectrum sources (e.g., radio galaxies and dif¬ 
fuse synchrotron sources) in particular will be strongly enhanced 
by MSSS. It is clear from Fig.[^ however, that a larger LOFAR 
array should be processed in order for the MSSS angular resolu¬ 
tion to be competitive with the other existing surveys. 


6.1. Long baselines 


A uniquely powerful aspect of LOFAR’s view of the low fre¬ 
quency sky is the sub-arcsecond resolution afforded by base- 
lines to and between international stations (see, e.g., [Varenius 
et al.||2015] ). International stations are always included in the 
LBA observations, but not in the HBA observations for the rea¬ 
sons described in § |2.2[ A long-baseline working group is plan¬ 
ning to use the same MSSS data that will generate the initial 
low-resolution images and ca talog as part of a sur vey for long- 
baseline calibrators (see also |Mold6n et al.||2Q15 ), and to pro¬ 
vide initial images for the brightest sources. Test processing 
rounds have demonstrated that higher angular resolution imag¬ 


ing is feasible despite the sparse uv coverage. HBA images at 5- 
30"resolution (using the outer Dutch remote stations) have been 
successfully produced. Efforts toward a large-scale reprocessing 
of the MSSS data to produce a high-resolution catalog are in an 
early phase at the time of writing. 


An example of higher resolution imaging of MSSS data 
is presented in Fig. The MSSS-HBA observations of held 
H063-r39 were reprocessed using an automated self-calibration 
loop that has been developed as part of an overall upgrade to 
the TIP (see § |3.4[ ). This scheme performs a progressively higher 
angular resolution sequence of source finding, phase calibration 
and imaging. Here, we have performed 8 iterations to achieve 
a final image resolution of 42" after averaging over 5 bands 
(cf. the initial resolution of MSSS images, 2.5' from baselines 
< 2kA; see § 3.4. 1| ). While the data provide the capability for 


still higher angular resolution, our aim here is to demonstrate 
image quality comparable to existing higher frequency radio sur¬ 
veys. We have selected 3 sources from the held of view that 
are resolved in the high resolution MSSS image: 3C 111 and 
two double-lobed radio galaxies (4C -r39.13 and 4C -r37.10). 
Through comparison with the NVSS images of the same objects 
it is clear that the morphology is properly reproduced in the re¬ 
processed MSSS data, at an order of magnitude lower frequency. 
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Frequency [MHz] 



Fig. 17. Comparisons between MSSS sensitivities (left) and resolutions (right) with those of other existing radio surveys as sum¬ 
marized in TableIn the lefthand panel, dashed lines indicate representative spectral indices of a = -0.7 and a = -\ A. The solid 
black lines illustrate the frequency dependence of the sensitivity in the 8 bands in each of the LB A and HBA survey segments, 
while the black stars show the frequency-averaged sensitivity demonstrated in §[^ In the righthand panel, the downward-pointing 
arrows indicate that the angular resolution of the initial MSSS catalog is limited with respect to the capabilities of the visibility data. 
Processing the full array will improve the survey performance. 





4h01nn 4h01m 4h01nn 4h01nn 


RA (J2000) 


Fig. 18. Images of three sources, one per row, identified in the field of H063-F39. The sources are 3C 111 (top), 4C -f39.13 (middle), 
and 4C -f37.10 (bottom). In each row the left p anel displays the VLSSr image (resolution 75"), the next two columns display MSSS 
images at 2.5' and 42" resolution (see § |6.1| ), and the right panel displays the NVSS image (resolution 45"). The beamsize for 
each column is plotted in the top row. Contour levels begin at 0.6 Jy beam“^ for VLSS; 1, 0.2 and 0.2 Jybeam"^ for low-resolution 
MSSS; 0.2, 0.15 and 0.1 Jybeam"^ for high-resolution MSSS; and 0.02 Jybeam"^ for NVSS. Contours increase by multiples of 2 
in all frames. 
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6.2. Transients and variable sources 


There are multiple aspects of the MSSS survey that are useful for 
conducting transient searches. First, the 9 (2) snapshots in LB A 
(HBA) mode allow a limited variability and transient capabil¬ 
ity. Furthermore, in LBA mode a single subband beam is always 
placed on the north celestial pole (NCP), enabling transient mon¬ 
itoring on longer timescales (albeit with reduced sensitivity). A 
full description of the NCP transient monitoring program, and 
the first low-frequency transient identified therein, is presented 
in a companion paper ( |Stewart|2015 1. 

MSSS is also expected to be useful to detect pulsars that 
are highly scattered to a degree that causes the pulse profile to 
be indistinguishable in beamformed observations (see Stappers| 
|et al.|20ir[ for a description of the beamformed pulsar observing 
method). In imaging observations, the total fiux density from the 
pulsar will still be visible, permitting study of the low-frequency 
spectral behavior of such objects. 



Faraday depth (rad m“^) 


6.3. Magnetism 


One of the fundamental scientific themes being pursued using 
LOFAR is the study of cosmic magnetism ( [Beck et al.|2013| ). A 
magnetism working group is p lanning to perform Faraday ro ta- 
tion measure (RM) Synthesis ( [Brentjens & de Bruyn||2005| ) on 
MSSS data in order to search for polarized sources, primarily 
Galactic (pulsars and diffuse foreground emission; see |Jelic et al. 
2Q14| ), but with h opes of detecting pola rized active galactic nu¬ 
clei (AGNs) (e.g., |Mulcahy et al.|2014| ) and giant radio galaxies 
(GRGs) as well. 

A polarization survey based on MSSS data will be uniquely 
powerful at these frequencies; with the full achievable angu¬ 
lar resolution and sensitivity over the entire survey area, it will 
help greatly in furthering our understanding of polarization at 
low frequencies. The HBA component of MSSS provides an 
excellent Faraday resolution of 1.3 rad m“^ and a maximum 
Faraday depth of approximately 330 rad m“^. The rotation mea¬ 
sure spread function (RMSF), which displays the response to 
simple polarized sources using all 8 HBA bands together from 
MSSS, can be seen in Fig. Initial tests of the polariza¬ 
tion recoverable through MSSS have already been performed. 
Polarized Galactic foreground emission from the “Fan region” 
has been detected, matching structures seen in previous obser¬ 
vations with WSRT (e.g., |Iacobelli et aL]2013| ). In addition to 
this, the highly polarized pulsar PSR J0218-F4232 with a known 


Faraday depth of -61 radm (Navarro et al. 1995) was de 


tected at the correct Far aday depth after correction for iono¬ 
spheric Faraday rotation ( [Sotomayor-Beltran et al.|[2QT^ . This 
demonstrates that an accurate polarization survey of the Galactic 
foreground and extragalactic background sources is feasible with 
MSSS-HBA. This effort will be presented in a future paper. 


6.4. Individual extragalactic and Galactic targets 

The key unique aspect of MSSS is its intrinsic broadband nature 
that enables study of the spectral properties in various classes 
of sources. The primary sources of interest include galaxy clus¬ 
ters, star-forming galaxies, AGNs, and supemovae, among oth¬ 
ers. Several working groups have been formed to investigate pre¬ 
liminary samples of targets and determine their low frequency 
spectra. Particular scientific questions to be addressed include 
the spectral properties of galaxy cluster halos and the search 
for previously unknown diffuse emission; the spectral behavior 
of nearby star-forming galaxies and the prevalence of spectral 


Fig. 19. Rotation measure spread function provided by MSSS, 
considering the combination of all 8 HBA bands. The width of 
the main lobe is 1.3 rad m“^, excellent for recovery of small 
RMs and thus weak magnetic fields. 


turnovers at low frequencies; characterization of GRGs and the 
search for previously unknown objects; and the search for new 
supernova remnants and Pulsar Wind Nebulae (PWNe). 


7. Summary and outlook 

We have presented the motivation and setup of the 
Multifrequency Snapshot Sky Survey (MSSS), a ground¬ 
breaking new radio continuum survey being performed with 
LOFAR. The primary design goal is a moderate angular resolu¬ 
tion (~ 2'), moderate depth (~ lOmJy beam“^), but intrinsically 
broadband survey to populate the initial LOFAR source database 
that will be used for calibration purposes in forthcoming deeper 
observations. The survey is also fertile ground for new scientific 
research into the source population of the low-frequency sky, 
and in particular uniquely enables the study of the spectral 
characteristics of the brighter sources detected by MSSS. 

The quality of the forthcoming initial MSSS data release 
has been illustrated within the MSSS Verification Field (MVF), 
a representative 100 square degree area centered at (a,S) = 
(15^, 69°). We find that the survey area studied here is 90% com¬ 
plete above 100 mJy in the HBA with angular resolution 108", 
and above 550 mJy in the LBA with 166" resolution. Source po¬ 
sitions and fiux densities match well with existing radio surveys 
in these frequency ranges, albeit with possible small fiux scale 
offsets that require confirmation over a larger sample area. 

Subsequent MSSS data releases are anticipated, with the po¬ 
tential for improved multi-directional calibration (essential es¬ 
pecially for the LBA portion of MSSS) and higher angular reso¬ 
lution imaging. These will be accompanied by a number of pa¬ 
pers focusing on specific science projects enabled by the MSSS 
catalog (as outlined in §|^ such as low-frequency transients, po¬ 
larized sources, the spectral properties of nearby star-forming 
galaxies, and steep spectrum radio galaxies. 
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